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Abstract 

Background: Four of the five species of Telopea (Proteaceae) are distributed in a latitudinal replacement pattern on 
the south-eastern Australian mainland. In similar circumstances, a simple allopatric speciation model that identifies 
the origins of genetic isolation within temporal geographic separation is considered as the default model. However, 
secondary contact between differentiated lineages can result in similar distributional patterns to those arising from 
a process of parapatric speciation (where gene flow between lineages remains uninterrupted during differentiation). 
Our aim was to use the characteristic distributional patterns in Telopea to test whether it reflected the evolutionary 
models of allopatric or parapatric speciation. Using a combination of genetic evidence and environmental niche 
modelling, we focused on three main questions: do currently described geographic borders coincide with genetic 
and environmental boundaries; are there hybrid zones in areas of secondary contact between closely related 
species; did species distributions contract during the last glacial maximum resulting in distributional gaps even 
where overlap and hybridisation currently occur? 

Results: Total genomic DNA was extracted from 619 individuals sampled from 36 populations representing the four 
species. Seven nuclear microsatellites (nSSR) and six chloroplast microsatellites (cpSSR) were amplified across all 
populations. Genetic structure and the signature of admixture in overlap zones was described using the Bayesian 
clustering methods implemented in STUCTURE and NewHybrids respectively. Relationships between chlorotypes were 
reconstructed as a median-joining network. Environmental niche models were produced for all species using 
environmental parameters from both the present day and the last glacial maximum (LGM). 
The nSSR loci amplified a total of 1 54 alleles, while data for the cpSSR loci produced a network of six chlorotypes. 
STRUCTURE revealed an optimum number of five clusters corresponding to the four recognised species with the 
additional division of T. speciosissima into populations north and south of the Shoalhaven River valley. Unexpectedly, the 
northern disjunct population of T. oreades grouped with T. mongaensis and was identified as a hybrid swarm by the 
Bayesian assignment test implemented in NewHybrids. Present day and LGM environmental niche models differed 
dramatically, suggesting that distributions of all species had repeatedly expanded and contracted in response to 
Pleistocene climatic oscillations and confirming strongly marked historical distributional gaps among taxes. 

Conclusions: Genetic structure and bio-climatic modeling results are more consistent with a history of allopatric 
speciation followed by repeated episodes of secondary contact and localised hybridisation, rather than with parapatric 
speciation. This study on Telopea shows that the evidence for temporal exclusion of gene flow can be found even 
outside obvious geographical contexts, and that it is possible to make significant progress towards excluding parapatric 
speciation as a contributing evolutionary process. 

Keywords: Allopatry, Environmental niche models, Parapatry, Secondary hybridisation, Speciation, Proteaceae, Telopea, 
waratah 
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Background 

Models of speciation are generally categorised on the role 
and nature of geographic isolation during the origination of 
reproductive barriers between differentiating populations. 
Usually, allopatric speciation is considered as the default 
model since it specifies a simpler, more biologically plausible 
process than other models [1,2]. In an allopatric model of 
speciation, genetic isolation resulting from temporal geo- 
graphic separation can lead to the evolution of reproductive 
barriers and speciation. Consequently, biogeographic pat- 
terns are of particular interest when attempting to test allo- 
patric models [3]. In contrast, the parapatric and sympatric 
speciation models involve the differentiation of lineages des- 
pite the possibility for ongoing exchange of genetic material. 
In these models, selective barriers must have a strong im- 
pact in order to preserve between-lineage divergence [4,5]. 

The spatial arrangement of genetically divergent popula- 
tions is important in allopatric speciation, as it involves the 
geographic separation of a previously continuous metapo- 
pulation by geological or climatic causes, or a founder event 
resulting from long distance dispersal [6]. Lack of gene 
flow, drift and secondary adaptation to local environmental 
conditions among geographically separated populations can 
eventually result in morphologically and genetically distinct 
lineages. However, current isolating mechanisms are not 
necessarily the same as those that initiated speciation and, 
as a result, secondary contact and hybridisation can ensue 
[3]. Although the presence of hybrid fronts between closely 
related lineages can confirm allopatric speciation with sec- 
ondary contact, such fronts are difficult to differentiate 
from zones of primary parapatric speciation. In fact, any 
mechanism that can cause divergence in allopatry can the- 
oretically also act in parapatry, as long as the selective gra- 
dient acting on differentiation is sufficiently strong to 
counterbalance continuing gene flow [1]. 

Here, we aim to identify the main speciation mechanism 
involved in the diversification of the plant genus Telopea 
(the waratahs; Proteaceae). Telopea is a genus comprising 
five species of long-lived, lignotuberous, sclerophyllous, 
bird-pollinated, wind-dispersed shrubs and small trees dis- 
tributed in mesic environments of south eastern Australia, 
including Tasmania. They live in eucalypt-dominated for- 
ests and woodlands, growing on acidic, nutrient-poor, well 
drained soils derived from sandstone, granite, chert or 
quartzite from sea level in the Sydney region to 1200 m alti- 
tude in subalpine habitats in Tasmania. Weston & Crisp [7] 
reconstructed Telopea phylogeny in a cladistic analysis of 
morphological characters and the species limits of the four 
mainland Australian species were tested morphometrically 
by Crisp & Weston [8]. 

Phylogenetically, the four recognised mainland taxa are 
grouped into two morphologically distinct clades comprising 
T, speciosissima and T. aspera to the north, and T. oreades 
and T. mongaensis to the south. Their distribution follows a 



latitudinal pattern ranging from southern Victoria to north- 
ern New South Wales (Figure 1). The species are mostly 
separated by large distributional gaps. The exceptions are a 
site at Monga National Park (NP), where a geographically 
disjunct population of T. oreades is found in sympatry with 
T. mongaensis, and an area also south of the Sydney Basin 
where a potential distributional overlap between the latter 
species and T. speciosissima occurs (recent fieldwork has 
located a single morphologically intermediate individual). 

Although terminally-winged, the seeds of Telopea are not 
easily dispersed and tend to germinate in proximity to the 
maternal parent [9]. As a result, the latitudinal distribution 
of the genus suggests that interruptions of connectivity be- 
tween once continuously distributed lineages could have 
contributed to allopatric speciation with the development of 
morphological and adaptive differentiation likely to have 
been secondary events. In a previous landscape genetic 
study on T. speciosissima, we suggested that substrate- 
mediated allopatry and temperature-dependent reproductive 
barriers were likely to affect temporal changes in connectiv- 
ity among populations [10]. While genetic differentiation in 
the south corresponded to a distributional gap caused by 
changes in edaphic conditions, altitudinal constraints in 
gene flow were the consequence of a strong association be- 
tween temperature and flowering time. 

In south-eastern Australia, the glacial cycles of the Qua- 
ternary shaped the distribution of plants and plant commu- 
nities by increasing aridity, fire frequency and fire intensity 
[11]. The limited fossil record suggests that dynamic succes- 
sional processes affected local vegetation types, and that 
localised moist refugia along the Great Dividing Range pro- 
vided some protection against localised extinctions [12,13]. 
Plant distributional changes during the Quaternary have 
been particularly well documented in northern latitudes 
[14,15] and more recently in eastern Australia [16]. Similar 
temporal changes can also be expected for Telopea. The ex- 
tent to which a species' distribution can respond to changing 
climatic conditions is affected by niche breadth, whilst the 
ability of a species to establish in new climatically suitable 
areas is dependant upon numerous ecological factors, in- 
cluding dispersal ability and competitive advantage. The use 
of environmental niche modeling (ENM) to investigate his- 
torical environmental suitability in combination with mo- 
lecular data can enhance our understanding of temporal 
changes in population dynamics [17]. Here we integrate 
genetic and bio-climatic methods to investigate speciation 
across all mainland species of Telopea and test the alterna- 
tive scenarios of parapatric or allopatric speciation. In par- 
ticular we focus on three main questions: 

1. Do currently described geographic borders coincide 
with genetic and environmental boundaries? 
Distributional and morphological data mostly 
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Figure 1 Distribution map for the four Telopea species investigated in this study. Map showing all verified distributional records (used as 
training data for the Maxent models) for the four study species (purple T. aspera; green T. speciosissima; red T. mongaensis; blue T. oreades). The 
major geographic points of interest discussed in the text are presented. 



support this [8] . Here we aim to validate between- 
species distinctiveness by combining present-day 
genetic (nSSRs, cpSSRs) and bio-climatic (ENM) 
analyses. 

2. Are there hybrid zones in areas of secondary contact 
between closely related species? Here we use targeted 
population genetic data in overlap zones to 
investigate hybridisation between species. A previous 
study on T. speciosissima suggested that admixture 
between previously differentiated populations can 
take place as a consequence of the post-LGM re- 
establishment of gene flow caused by temperature- 
mediated phenological shifts along the altitudinal 
gradients [10]. 



3. Did species distribution contract during the last 
glacial maximum, resulting in distributional gaps 
even where overlap and hybridisation currently 
occur? We use ENM to predict historical ranges for 
the four species and gather evidence for 
distributional discontinuities. The existence of such 
patterns would provide support for secondary rather 
than primary contact zones and consequently favour 
a scenario of allopatric over parapatric speciation [3] . 

Methods 

Study species 

Five species of Telopea are distributed from the Gibraltar 
Range in south-eastern mainland Australia (29° 28' S) to 
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southern Tasmania (43° 33' S). Of the four mainland species, 
T. aspera in the north is the only one completely separated 
from its nearest congeners by about 400 km, while the other 
three have potential contact points at the extremes of their 
distributions (Figure 1). All species have naturally fragmen- 
ted distributions, of which the most intriguing is that of 
T. oreades, the most northern population of which overlaps 
with the most southerly population of T. mongaensis but is 
located over 160 km from the next most northerly popula- 
tion of T. oreades [8]. Telopea mongaensis also shows a sig- 
nificant internal gap of over 45 km imposed by the 
Shoalhaven River valley, with its most north-easterly popula- 
tions occurring within 20 km of the nearest populations of 
T. speciosissima. Telopea speciosissima is the most wide- 
spread species with populations south of the Shoalhaven 
River valley being isolated from those on the northern side 
by a over 30 km [10]. 

According to the morphology-based cladistic analysis of 
Proteaceae subtribe Embothriinae published by Weston & 
Crisp [7], a clade of T. speciosissima and T. aspera is sister to 
a clade including T. truncata (Tasmania), T. oreades and 
T. mongaensis. The speciosissima-aspera clade is characterised 
by three macromorphological synapomorphies: presence of 
marginal leaf teeth (versus entire or lobed margins), and 
enlarged, bright red involucral bracts (versus small and green, 
to dusky pink involucral bracts). The truncata-oreades-mon- 
gaensis clade is characterised by basitonic inflorescences 
(flowers open from the tip of the inflorescence towards the 
base, in contrast to the acrotonic inflorescences of the specio- 
sissima-aspera clade and outgroups, in which flowers open 
from the base towards the tip). Telopea oreades and T. mon- 
gaensis form a clade characterised by markedly bicolorous 
perianths in which the adaxial surfaces are a much brighter 
shade of red than the abaxial surfaces, in contrast to the con- 
colorous perianths in other Telopea species. 

Sampling, DNA extractions and PCRs 

Sampling was aimed at obtaining an account of the genetic 
diversity across the entire distribution of mainland species 
(Table 1). Leaf material was collected from a total of 619 
individuals, and total genomic DNA was extracted using 
DNeasy® 96 plant kits (QIAGEN®, Hilden, Germany). 

Seven nuclear microsatellites (nSSR) specifically developed 
for T. speciosissima (TS03bgt, TS04bgt, TS12bgt, TS13bgt, 
TS18bgt, TS23bgt and TS27bgt) were amplified across all 
individuals using the PCR protocols of Porter et al. [18]. 
Genotyping results were checked with Microchecker [19] 
for evidence of scoring errors due to stuttering, large allele 
dropout and null alleles with ranked-based and binomial- 
based estimates showing no evidence of scoring error. 

Chloroplast microsatellites (cpSSR) were amplified from 
three individuals per population, using ccmplO [20], psbA 
[21], trnG (UCC)exl - atpA [22], Lomcpl, Lomcp2 and 
Lomcp3 (M Milner pers. com.). Reaction conditions were 



as described in Ebert & Peakall [22] using Schuelke's [23] 
method of fluorescently-labelled M13 forward primers. 
Amplification products were visualized using a 3130x1 
Genetic Analyzer (Applied Biosystems) and scored for size 
using GeneMapper 3.7 (Applied Biosystems). 

To test genotyping accuracy approximately 20% of all PCR 
reactions were repeated for each primer pair across random 
individuals. Less than 5% of those PCR repeats identified 
errors that needed confirmation through replication of PCR 



Table 1 Study sites 





Species 


Population 


N 


Latitude 


Longitude 


1 


T. aspera 


Washpool NP 


20 


29°29.11'S 


152°19.39'E 


2 


T. aspera 


Gibraltar Range 


22 


29°30.8'S 


152°21.75'E 


3 


T. aspera 


Mulligans Rd 


20 


29°31.73'S 


152°20.19'E 


4 


T. aspera 


Anvil Rock 


22 


29°33.1'S 


152°18.88'E 


5 


T. speciosissima 


Kulnura 


21 


33°12.92'S 


151°1 1.92'E 


6 


T. speciosissima 


Newnes Forest 


20 


33°23.67'S 


150°12.78'E 


7 


T. speciosissima 


Mountain Lagoon 


20 


33°26.82'S 


150°38.55'E 


8 


T. speciosissima 


Bell Line Rd 


20 


33°29.93'S 


150°16.08'E 


9 


T. speciosissima 


Kurrajong Heights 


10 


33°31.13'S 


1 50°37.22'E 


10 


T. speciosissima 


Patonga 


19 


33°32.23'S 


151°1 7.05'E 


11 


T. speciosissima 


West Head 


15 


33°36.57'S 


15H6.52'E 


12 


T. speciosissima 


Duffys Forest 


22 


33°39.68'S 


151°11.5'E 


13 


T. speciosissima 


Watagan SF 


21 


33°4.37'S 


151°20.05'E 


14 


T. speciosissima 


Kings Tableland 


21 


33°45.97'S 


150°22.93'E 


15 


T. speciosissima 


Mt. Alexandra 


20 


34°26.58'S 


150°27.23'E 


16 


T. speciosissima 


Carrington Falls 


19 


34°37.53'S 


150°39.32'E 


17 


T. speciosissima 


Brogers Creek 


8 


34°42.00'S 


150°40.60'E 


18 


T. speciosissima 


Waterfall Flat 


15 


34°8.83'S 


151°0.33'E 


19 


T. speciosissima 


Bottle Forest 


11 


34°8.92'S 


151°4.62'E 


20 


T. speciosissima 


Curra Moors 


10 


34°8.92'S 


151°4.62'E 


21 


T. speciosissima 


Ulladulla 


21 


35°22.15'S 


150°28.4'E 


22 


T. speciosissima 


Turpentine Range 


21 


35°3.72'S 


150°25.42'E 


23 


T. speciosissima 


Jervis Bay 


20 


357.93'S 


150°41.1'E 


24 


T. mongaensis 


Gunrock Creek 


10 


34°37.75 , S 


150°24.72'E 


25 


T. mongaensis 


Budawangs 


20 


35°24.74'S 


150°1.16'E 


26 


T. mongaensis 


Dasyurus PA 


20 


35°33.65'S 


149°55.32'E 


27 


T. mongaensis 


River Forest Rd 


11 


35°37.05'S 


149°54.7'E 


28 


T. mongaensis 


Monga NP 


10 


35°37.93'S 


149°54.73'E 


29 


T. oreades 


Monga NP 


29 


35°37.93'S 


149°54.73'E 


30 


T. oreades 


Waratah Creek 


20 


37°0.31'S 


149°22.75'E 


31 


T. oreades 


Errinundra NP 


20 


37°16.51'S 


148°53.23'E 


32 


T. oreades 


Bonang Rd 


19 


37°24.98'S 


148°36.18'E 


33 


T. oreades 


Cooaggalah Rd 


10 


37°26.06'S 


149°19.99'E 


34 


T. oreades 


Combeinbar Trail 


20 


37°28.29'S 


149°5.91'E 


35 


T. oreades 


Lind NP 


5 


37°34.39'S 


148°58.02'E 


36 


T. oreades 


Mt Bemm 


7 


37°38.82'S 


148°58.71'E 



Name, sample size (N) and location of each Telopea population included in 
this study. 
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and genotyping. Finally, to verify that the amplified loci in- 
deed contained the expected nSSR repeats across all Telopea 
species and to verify the accuracy of allelic sizes, representa- 
tives of species/locus combinations were sequenced [Gen- 
Bank accession no. JF931649-JF931672]. Unique cpSSR 
allele sizes were also sequenced across species to ascertain 
that the variation was due to mononucleotide repeat regions 
[GenBank accession no. JQ778988-JQ779019]. 

Defining genetic boundaries 

Species-level averages of nSSR allelic diversity (A), expected 
(H e ) and observed (H D ) heterozygosities as well as number 
of private alleles (A p ) were calculated using GDA 1.1 [24]. 
In order to avoid bias caused by uneven sampling [25], a 
standardized estimate of allelic richness (R63) independent 
of sample size [26] was calculated using the program 
FSTAT 2.9.3 [27]. Pariwise F st [28] values were also 
obtained through FSTAT and used to provide within and 
between species averages. 

We used the Bayesian clustering method described by 
Pritchard et al. [29] and implemented in STUCTURE 2.3.3 
on the nSSR data to identify genetically differentiated groups 
of individuals in the absence of preliminary information on 
group boundaries. The model assumes the existence of K 
clusters (the real number being unknown) and uses the al- 
lelic frequencies at each locus to assign individuals to these 
clusters through a Markov chain Monte Carlo (MCMC) 
probabilistic approach in which individuals are assigned to 
clusters so as to maximize Hardy- Weinberg equilibrium 
within populations. All analyses were based on 10 independ- 
ent runs for each K value, with each individual run being 
based on 2.5x10 s MCMC iterations following a burn-in 
period of 5x10 iterations without prior information on the 
taxonomy or the locality of origin of the individuals sampled. 
The admixture frequency model was run under the assump- 
tion of correlated allele frequencies to improve clustering of 
related lineages and identify possible hybridisation patterns 
[30]. The optimal number of clusters was verified using the 
AK statistical approach suggested by Evanno et al. [31]. 

After a preliminary test aimed at finding a suitable range 
for K and the optimal burn-in period, we tested K from 1 to 
20 for the full dataset including the four species under study. 
We also tested K from 1 to 12 for the dataset including only 
T. oreades and T. mongaensis to specifically investigate as- 
signment patterns for the only sympatric individuals located 
at the Monga NP site. 

Analysis of molecular variance (AMOVA; [32]) was then 
used to quantify variance components and the significance 
of the genetic subdivisions identified by the Bayesian tests, as 
well as a range of other relevant groupings. Permutation tests 
for two-level AMOVAs were implemented to test whether 
levels of differentiation were significantly greater than zero. 
AMOVA results were obtained using GenAlEx6.3 [33]. 



Assuming that cpDNA does not recombine and can be 
treated as a single locus, we combined the sizes of all six 
cpSSR fragments to give the one chloroplast haplotype 
(chlorotype) per individual. A median-joining network [34] 
was constructed using NETWORK v4.5.1.6 [35], a program 
that constructs networks based on size differences among 
haplotypes. 

Detecting hybridisation 

Since STUCTURE (and the other analytical approaches 
used) is not explicitly designed to assign individuals to 
admixed classes, we used the Bayesian model implemented 
in NewHybrids [36] to detect admixed individuals resulting 
from the interbreeding between distinct species. The ori- 
ginal tests and simulations by Anderson & Thompson [36] 
show that admixture can be detected without the need for 
diagnostic alleles. Although a high number of informative 
loci produces much better posterior probabilities (PP) for 
assigning to hybrid categories, when F st > 0.2 (as is the case 
in this study; Table 2) a smaller number of loci is sufficient 
for detecting hybrids [37]. 

Our objective was not to identify and quantify specific ad- 
mixture categories but rather to ascertain whether differ- 
entiated genomes, particularly those of T. oreades and 
T. mongaensis, have been mixed where the species overlap. 
To this end, we conducted a series of preliminary analyses 
that culminated in a run using all 11 populations from the 
two species. Jeffreys priors were used with a burn-in of 
1x10 s sweeps followed by 3.5x10 s sweeps. Posterior prob- 
ability of assignment as pure, Fls, F2s, and backcrosses 
were initially measured, and proportions of admixture vs. 
pure were then calculated from these results. An individual 
was assigned as a pure or hybrid individual when its poster- 
ior probability of belonging to that genotypic class exceeded 
PP = 90% [38,39]. 

Factorial correspondence analysis (FCA) of diploid gen- 
etic data [40] was carried out using GENETIX 4.05 [41] to 
produce comparative graphical representations of genotypic 
relationships between all T. oreades and T. mongaensis indi- 
viduals sampled. 

NewHybrids was also used to analyse the T. speciosissima 
and T. mongaensis populations that were geographically 
closest to the only known morphological intermediate be- 
tween these two species. This was to serve as a test for the 
analytical approach and to investigate the possible presence 
of introgression. 

Current and historical environmental niche modelling 

We used the machine-learning automated statistical algo- 
rithm Maxent 3.3.3 (Maximum entropy modeling applica- 
tion; [42]) to identify the multivariate correlations between 
the 19 available -920 m x 770 m pixel (30 arc-second) 
WorldClim [43] climate axes (environmental data) and cli- 
mate conditions present at geographic localities matching 
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Table 2 Diversity values for the taxa investigated 





Pops 


N 


A 


Ad 




Mean H 

i v i cd ii i Iq 


Mean H 

1 V 1 CO III Iq 


F S f wtn 


F st btw 


Mpan Ci 

IVlCul 1 


cpSSR 


T. aspera 


4 


84 


9.3 


18 


9 


0.554 


0.626 


0.046 


0.265 


0.95 


1,2 


T. mongaensis 


5 


70 


6.3 


7 


6.3 


0.38 


0.492 


0.192 


0.318 


0.96 


4,6 


T. mongaensis plus T. oreades at Monga NP 


6 


98 


6.7 


7 


6.9 


0.415 


0.542 


0.22 


0.303 


0.95 


4,6 


T. oreades 


8 


129 


7.3 


12 


6.8 


0.482 


0.619 


0.238 


0.252 


0.84 


4 


T. oreades minus T. oreades at Monga NP 


7 


98 


6.5 


11 


6.3 


0.472 


0.564 


0.201 


0.276 


0.95 


4 


T. speciosissima 


19 


334 


17.2 


54 


12.7 


0.686 


0.794 


0.112 


0.222 


0.95 


3,4,5,6 


T. speciosissima south of Shoalhaven only 


3 


62 


8.5 


6 


7.6 


0.648 


0.714 


0.093 


0.243 


0.90 


6 



The rows report species (or groups pf population) specific values. Pops: number of populations; N: number of individuals; A: allelic diversity; A p : number of private 
alleles for the species; R 63 : unbiased allelic diversity; Mean H e and Ho: mean expected and observed heterozygosities; F st wtn: average F st value among 
populations within a species; F st btw: average F st value relative to other species; Mean Q: mean individual-level assignment value to a species (based on the 
Bayesian clustering approach implemented in STRUCTURE, at K = 6); cpSSR: identity of chloroplast haplotypes amplified. 



the known native distributions of each of the Telopea spe- 
cies. We employed this environmental niche modeling 
(ENM) technique to determine whether present day habitat 
characteristics correspond to genetic divergence between 
taxa. Maxent was chosen as it has been shown to outper- 
form other modeling methods when generating predictions 
of species' ranges [44]. 

Occurrence records were compiled from all Australian 
Herbaria (voucher - high taxonomic integrity) and taxo- 
nomically corrected records from the (non-voucher - high 
geographic integrity) databases: YETI plots [45]; and 
ATLAS incidental observations [46]. We ran multiple 
model sets for Telopea as a whole and for each Telopea 
species to determine consistent geographic extents between 
models and consistently dominant climatic factors [47]. 
Model run 1) only included locality records with accuracies 
better than 100 m (Table 3). These high precision sites were 
also of high accuracy (giving ±0.1 climatic pixel accuracy) 
and those sites that fell within 100 m of the climatic pixel 
boundary were excluded. This gave a single climatic pixel 
that best represents the sites sampled. This model repre- 
sents a minimum estimate of the realised climatic space 
(climatic conditions in which the species is specifically 
recorded) of each species but could result in over-fitting of 
the models [47]. Model 2) included Herbaria records of 
geographic accuracies better than 1 km, (±1 climatic pixel 
accuracy or a potential of 9 climatic pixel values per site); 
and model 3) the above plus non-voucher records better 
than 1 km accuracy (Table 3). Whilst the climatic data sub- 
mitted to the model is less precise, substantially more 
points are submitted providing more robust correlations. 
Model runs 2) and particularly 3) provide better estimates 
of the fundamental climatic space (the climatic conditions 
under which the species is capable of reproducing and per- 
sisting) and minimise the effect of false omissions from the 
training localities which would usually result in models 
underestimating suitable climatic ranges [47]. The three 
sets of model runs were compared for between-model con- 
sistencies. The larger numbers of localities submitted to the 



training data better capture the realised climate space and 
avoid having (climatic) variable numbers approaching the 
number of sample locations which leads to over-fitting of 
the models [47]. 

Model outputs were projected onto to the eastern states 
of Australia using the same 19 WorldClim variables for 
0 ka (pre-industrial) and 21 ka (LGM) time frames. Models 
were run in the default settings but with 10 replicates under 
the cross-validation option, and with response curves and 
jackknife settings enabled. However, this resulted in over- 
fitted models (as shown geographically by over fragmenta- 
tion and underestimation of suitable climatic habitat within 
expected available soil/plant community types) so the mod- 
els were rerun with hinge features disabled. 

The models were run initially with all 19 variables to deter- 
mine the relative contributions of each environmental vari- 
able to the Maxent model. Then, to reduce the effects of 
overfitting [47], only variables that were not auto-correlated 
and contributed more than 5% to the Maxent model were 
included in a subsequent model run. 

Telopea species are highly constrained by soil type, and 
non-climatic abiotic factors like this can constrain realised 
distributions to portions of the fundamental climatic space. 
If obviously dominant factors of this nature are excluded 
from the assessment of suitable geographic extents to sub- 
mit to the training of the models, then this has the effect of 
submitting more false omissions to the model, producing 
an underestimation of suitable habitat [47]. To better deter- 
mine the available realised habitat and estimate the (true) 
pseudo-absences within the fundamental space, the geo- 
graphic extents of the training data were clipped to suitable 
soil types. To make soil type allocation relevant to current 
and historic time frames, surface geological units were used 
rather than soil type and were buffered to 1 km to capture 
likely surface soil types. Soil types for each species were 
assessed by overlaying occurrence records with soils map- 
ping and resulted in the identification of mostly sediment- 
ary siliciclastic, quartz-rich arenite to rudite or igneous 
felsic intrusive lithologies. 
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Table 3 Site details for the ENM models 





Model 1: 100 m Herbarium voucher 


Model 2: 1 km Herbarium voucher 


Model 3: 1 km Herbarium and non-voucher 


T. aspera 


4 


14 


113 


T. mongaensis 


4 


24 


150 


T. oreades 


7 


37 


259 


T. speciosissima 


19 


85 


647 



Numbers of sites analysed for each of the four Telopea species per ENM run type. 



In order to maintain the integrity of the pseudo-absence 
reference that is integral to the software and to reduce the 
effects of over fitting [47], model runs were repeated for the 
19 training environmental variables from 0 ka (pre-industrial 
time) and the geographic extents of training data were con- 
strained to within 200 km of: 1) the genus point locations; 2) 
the species point locations; 3) the species point locations and 
to suitable soil type. These three sets of model runs represent 
increasing removal of false omissions from the training data. 
Each of these model runs was projected onto the following 
two geographic extents: a) the eastern Australian states; and 
b) the eastern Australian states constrained by available suit- 
able soil type. The latter output represented a reduction in 
the false commissions in the projection. 

To determine the extent of overlap between species pairs, 
the outputs of the last model run (3b) were intersected for 
each species pair in ArcGIS 9.2 to produce a map layer of 
the probability of the two species to overlap (Table 4). 

Results 

Defining genetic boundaries among lineages 

The seven SSR loci amplified a total of 154 alleles across the 
36 populations representing four recognised Telopea species. 
Species-level allelic richness ranged from Rg 3 = 6.3 to 17.2, 
and heterozygosity measures from H e = 0.492 to 0.794 and 
from Hq = 0.380 to 0.686 (Table 2). Telopea speciosissima 
was consistently the most, and T. mongaensis consistently 
the least diverse species. When analysed on their own, the 



three most southerly populations of T. speciosissima still 
contained high levels of diversity, while shifting the putative 
T. oreades population at Monga NP to T. mongaensis made 
little difference to their respective diversity measures 
(Table 2). 

Average pairwise F st values suggest strongest genetic dif- 
ferentiation within T. oreades (particularly when the Monga 
NP population is included) and weakest within T. aspera. 
Telopea mongaensis was the most genetically differentiated 
from all other taxa (Table 2). 

Using the seven nSSR loci, STRUCTURE'S Bayesian clus- 
tering produced the most substantial increases in L«P(D) 
values at K = 5, with the AA" statistic preferentially supporting 
K = 5 (AKs = 358.7) as the most robust division. Here we 
visualise data from the next most robust division (AK 6 
= 198.5) as it better represents previously-described sub- 
structure within T. speciosissima [10]. Assignment values of 
Q > 95% were usually obtained at the species level (Table 2). 
However, two major taxonomic inconsistencies emerged: a) 
the assignment of the T oreades population at Monga NP to 
T. mongaensis; b) the separation of the southern populations 
of T. speciosissima from the rest of the species. 

AMOVA supported significant segregation among Telopea 
species. High levels of genetic variance were partitioned to 
between-species differentiation (42.7% of variation; P < 0.001), 
although similar levels of genetic differentiation were parti- 
tioned to within-species (43.6% of variation; P < 0.001). A total 
of six distinct chlorotypes were obtained across the four 



Table 4 Modelled current (pre-industrial) and LGM distribution overlap among four Telopea species 



Species 






T. aspera 




T. speciosissima 




T. mongaensis 






P 


Area (ha) (p>=0.05) 


P 


Area (ha) (p > =0.05) 


P 


Area (ha) (p > =0.05) 


T. speciosissima 


LGM 


0.142 


141 












Current 


0.247 


2400 












% Change 


74 


1603 










T. mongaensis 


LGM 


0.002 


ND 


0.361 


213299 








Current 


0.004 


ND 


0.616 


498960 








% Change 


100 


ND 


71 


134 






T. oreades 


LGM 


0.033 


ND 


0.391 


6782 


0.535 


142955 




Current 


0.039 


ND 


0.452 


120901 


0.202 


22041 




% Change 


18 


ND 


16 


1683 


-62 


-85 



p represents the maximum modeled likelihood of geographic distribution overlap between species pairs within any one pixel. The area (in hectares) represents 
overlap areas of maximum modeled likelihood of overlap exceeding 5% for both species in the comparison.% Change is the change in species overlap extent 
from LGM to current. ND is no area detected for p > = 0.001 by the 30 arc second pixel models. 
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mainland species and none of the sampled populations ampli- 
fied more than one chlorotype each (Table 2). A single, shared 
chlorotype (chlorotype 4) was amplified for all T. oreades 
(Figure 2g), all but one population of T. mongaensis 
(Figure 2e) as well as the Brogers Creek population of T. spe- 
ciosissima (Figure 2c). Multiple chlorotypes were obtained for 
T. speciosissima, with chlorotype 6 being found only among 
the southern populations as well as the Gunrock Creek popu- 
lation of T. mongaensis (Figure 2e). 

Detecting hybridisation 

Further STRUCTURE analyses were conducted to refine 
and clarify the interpretations of the assignment among the 
two southern taxa. Bayesian clustering across 199 indivi- 
duals from the two southern Telopea species, produced 
substantial increases in L«P(D) values at K = 2 and K = 3, 
with the AK statistic preferentially supporting K = 3 (AK 2 
= 417 .4; AK 3 = 1009.7). This result identified the southern- 



most population at Monga NP (including both putative 
T. oreades and T. mongaensis individuals) as differentiated 
from the two species. K = 2 differentiated between the two 
species but assigned the T. oreades population at Monga 
NP to T. mongaensis. The FCA plot provided further 
insight in the K = 3 result by placing the individuals at the 
Monga NP site (morphologically identified as either of the 
parental species or hybrids) as a third intermediate group 
suggesting either a differentiated lineage or a hybrid popu- 
lation (Figure 3a). 

The Bayesian assignment test implemented in NewHy- 
brids correctly identified the only morphologically intermedi- 
ate individual between T. speciosissima and T. mongaensis 
(collected at Brogers Creek) as a hybrid (PP = 92%) when no 
priors were given, and as an Fl hybrid (PP = 98%) if an Fl 
prior was used in the analysis. No further sign of admixture 
was detected between the two species, even among the 
other individuals at Brogers Creek (average assignment to 








Figure 2 Modelled bio-climatic distribution of four Telopea species (LGM and current) and measured population-level diversity 
(chloroplast and nuclear). Maxent modeled bio-climatic distribution for each Telopea species for the LGM (21 ka) and current (pre-industrial) 
time frames (showing only model run 3b being all records projected only onto soils suitable for Telopea), showing population-level nSSR and 
cpSSR (colour-coded haplotypes) data for each of the sampled populations. The nSSR pie charts partition mean individual-level assignment value 
to a species (based on the Bayesian clustering approach implemented in STRUCTURE, at K=6). The chloroplast network is shown in (a), and the 
model probability (likelihood occurrence of suitable habitat) is in a gray scale from 0 (absent) to1 (certain presence) in (b). a) T. aspera, LGM 
habitat and cpSSR; b) T. aspera, current habitat and nSSR; c) T. speciosissima, LGM habitat and cpSSR; d) T. speciosissima, current habitat and nSSR; 
e) T. mongaensis, LGM habitat and cpSSR; f) T. mongaensis, current habitat and nSSR; g) T. oreades, LGM habitat and cpSSR; h) T. oreades, current 
habitat and nSSR. 
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T. speciosissima for the other Brogers Creek individuals; 
PP = 98%). 

NewHybrids analysis of the T. mongaensis and T. oreades 
populations identified extensive hybridisation at the Monga 
NP site (Figure 3b), with 13% of individuals assigned to 
pure T. mongaensis (PP > 90%), and the remainder assigned 
to a potential hybrid category (PP < 90% of being pure) 
with 60% having PP > 90% of being either F2 or back-cross 
to T. mongaensis. None of the Monga NP individuals 
(identified based on morphology to either T. mongaensis or 
T. oreades) were assigned to pure T. oreades (only four 
individuals had a PP > 0% of being pure T. oreades). At the 
other sites, all individuals (with the exception of three 
T. mongaensis individuals at River Forest Rd, the site clos- 
est to the Monga NP site) were assigned with PP > 90% to 
either pure T. oreades or pure T. mongaensis. 

Current and historical environmental niche modelling 

The average test AUC for the 10 replicate runs indicates ex- 
cellent model fit for all species tested across all model sets 
indicating an excellent fit to estimated realised climatic 
spaces for Telopea and for each species tested [47]. For the 
model sets l-3b (trained on extents to within 200 km of the 
species point locations and to suitable soil type, and pro- 
jected onto only suitable soil types in the eastern Australian 
states) the average test AUC was: T. aspera 0.998 (std dev 
0.001); T. speciosissima 0.943 (std dev 0.018); T. mongaensis 
0.990 (std dev 0.006); and T. oreades 0.983 (std dev 0.011). 
However, this does not mean that the models sufficiendy 
estimated the fundamental climatic distributions [47]. The 
model run for the Herbarium records of T. aspera with 
100 m accuracy (with lower number of training localities) 
was the most fragmented, indicating overfitting of the model. 
Fragmentation of the remaining model sets within the 
expected distribution was minimal, particularly for the soils 
runs, indicating good representation of the fundamental 
climatic space for each species. Interestingly, the distribu- 
tion of soil types suitable to Telopea closely matches the 
modelled climatic distribution, consequently modelled cli- 
matic probabilities and areas did not differ appreciably be- 
tween model runs, at species or genus levels, with or 
without soils included. Telopea aspera responded most to 
higher precipitation in the warmest quarter. Telopea 
speciosissima responded to higher mean temperatures for 
the warmest and wettest quarters and to a higher precipita- 
tion in the driest quarter. Telopea mongaensis responded to 
higher temperatures in the wettest quarter and lower tem- 
peratures in the driest quarter (low seasonal temperature 
and precipitation variation). Telopea oreades responded to 
low maximum temperatures (for the warmest month and 
quarter). 

The modelled available bio-climatic habitats for T. aspera 
and T. speciosissima expanded from the LGM to present 
(Figure 2a-d). In contrast, the distribution of habitats 



suitable for T. mongaensis and T. oreades declined from the 
LGM to present (Figure 2e-h). As Telopea generally occurs 
in topographically variable habitats it is probable that there 
was greater in-situ preservation than the LGM modeling 
indicates, particularly in microhabitats too fine to be 
observed in the model pixel size of 30 arc seconds (approxi- 
mately 850 m). 

The LGM bio-climatic models for the various lineages 
support appreciable overlap only for T. speciosissima and 
T. mongaensis, and T. mongaensis and T. oreades. Due to 
the modelled current expansion of T. aspera and T. specio- 
sissima, overlap of distributions increased between all spe- 
cies pairs excluding between T. mongaensis and T. oreades 
(Table 4). Although some low-level potential overlaps are 
modelled between these two species, a 160 km wide distri- 
butional gap currently exists. Similarly, despite a large 
potential habitat overlap between T. speciosissima and 
T. oreades they are significantly separated geographically. 
Similarly T. aspera is currently too geographically disjunct 
(400 km) from all species for any interaction to occur. 

Discussion 

Concordance between geographic range and the 
distribution of genes and suitable habitat 

The species boundaries described for Telopea [7] follow a 
pattern of latitudinal differentiation along the Great Divid- 
ing Range (GDR). The GDR, running north-south along 
the length of eastern Australia, is a dominant topographic 
feature. Despite having relatively low elevation, it is suffi- 
ciently close to the coast to provide precipitation gradients. 
These gradients support a range of species and communi- 
ties that differ considerably from those found further west 
across the prevalently dry continent. The GDR supported 
refugial habitats during the climatic oscillations of the Qua- 
ternary, as revealed by regional species turnover [48,49] and 
deep between-population genetic divergences [50,51]. Based 
on palaeoecological evidence [12,13], phylogeographic diver- 
gences are interpreted as the genetic signatures of distri- 
butional contractions caused by climate-induced loss of 
suitable habitat. Habitat contractions were followed by re- 
expansions at the return of more favourable conditions [52- 
54]. These recurring Quaternary patterns have the potential 
to isolate populations and initiate allopatric speciation, even 
in the absence of strong intrinsic reproductive barriers 
[55,56]. 

The nSSR data provided fine-scale resolution and broad 
support for latitudinal species-level genetic divergences 
(Figure 2). The highest diversity was measured in the most 
widespread species (T! speciosissima) and, as expected in an 
allopatric speciation model [57], the more diverse species 
also showed signs of significant structure particularly in rela- 
tion to latitudinal landscape barriers and elevation. Current 
environment niche models are also in agreement with spe- 
cies being differentiated along the latitudinal gradient, 
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Figure 3 Testing hybridisation patterns at a distributional overlap zone (T. oreades and T. mongaensis at Monga NP). a) Factorial analysis 
comparing all individuals from all tested populations from T. oreades (white squares) and T. mongaensis (black triangles). The individuals from the 
Monga NP site morphologically assigned to either of the two species are represented in grey, b) NewHybrids results comparing the individuals 
(single bars) from the Monga NP site to the geographically closest populations of T. mongaensis (Dasyurus PA, River Forest Road) and T. oreades 
(Waratah Creek). T.m. DPA: T. mongaensis - Dasyurus PA; T.m. RFrd: T. mongaensis - River Forest Road; T.m. MNP: T. mongaensis - Monga NP 
(sympatric); T.o. MNP: T. oreades - Monga NP (sympatric); T.o. WC: T. oreades - Waratah Creek. Black represents% assignment to T. mongaensis 
genotype, white represents% assignment to T. oreades, grey represents% assignment to hybrid. 



suggesting relatively small bio-climatic overlaps between spe- 
cies (Table 4). Genetic and bio-climatic boundaries corres- 
pond to breaks in the GDR. T. aspera and T. speciosissima 



are separated by the Hunter River Corridor (an important 
floristic boundary; [58,59]). With the exception of the T. 
oreades population at Monga NP, T. mongaensis and T. 
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oreades are separated by an elevation decline of the GDR 
(previously reported as affecting genetic structure in other 
plants; [60]) and the fragmented nature of suitable edaphic 
conditions. 

Chloroplast SSRs could not fully resolve species boundar- 
ies (only six chlorotypes were discovered across the entire 
distribution of the four continental species), however the 
geographic distribution of the chlorotypes confirmed a lati- 
tudinal partitioning of genetic variation, even confirming 
within-species latitudinal divergences in T. speciosissima 
(Figure 2). Thus, genetic and bio-climatic analyses largely 
corroborate an allopatric speciation model where the lati- 
tudinal gradient is a major factor influencing lineage div- 
ergence. Nevertheless, present-day distributions include 
potential contact zones between species, and localised 
chlorotype-sharing suggests that between-species admixture 
is possible. If confirmed, hybridisation events could provide 
supporting evidence either for secondary contact between 
once-disjunct taxa or for parapatric differentiation along the 
borders separating sister taxa. Consequently, we explored in 
greater detail the possible occurrence of hybridisation at 
known contact zones. 

Supporting evidence for hybridisation at contact zones 

Artificial hybrids among Telopea species are common in 
the horticultural trade and all species are known to be 
potentially inter-fertile in cultivation [61]. A single and iso- 
lated individual that is morphologically intermediate be- 
tween T. speciosissima and T. mongaensis was confirmed as 
a hybrid (potentially an Fl) by the nSSR analyses. No 
T. mongaensis individuals were found in the vicinity, and 
the closest Telopea population was T. speciosissima at 
Brogers Creek. Interestingly, the Brogers Creek population 
shared chlorotype 4 with T. mongaensis despite bearing no 
morphological or nuclear signature of admixture. Further- 
more, while the nSSR data preferentially assigned this same 
population to the southern T. speciosissima group, the geo- 
graphically closest T. mongaensis population shared chloro- 
type 6 with the southern T. speciosissima, while being 
unequivocally assigned to T. mongaensis by nSSRs (Figure 2). 
These patterns suggest temporal variations in species overlap 
in response to changing climatic conditions, leaving genetic 
signatures similar to those observed in Tasmanian eucalypts 
[62,63]. 

Telopea oreades and T. mongaensis are sister taxa and 
the only qualitative morphological character that differenti- 
ates them is the absence of foliar sclereids in T. mongaensis. 
They also differ in average leaf dimensions, bract width and 
prominence of leaf venation but variation in these charac- 
ters overlaps between the species [8]. Consequently, mor- 
phological intermediates are difficult to characterise. These 
two species are generally geographically disjunct, but vari- 
ation in leaf morphology, including the presence of foliar 
sclereids in large-leaved plants, suggested that T. oreades 



and T. mongaensis may occur in sympatry at Monga NP 
(the southern distributional end of T. mongaensis). While 
both species share chlorotype 4 (Figure 2), which could be 
interpreted as a result of incomplete lineage sorting in view 
of limited chloroplast variation, the nSSR-based Bayesian 
analyses implemented in NewHybrids presented an unam- 
biguous admixture scenario at Monga NP (Figure 3) sug- 
gesting repeated backcrossing to T. mongaensis but no 
evidence of the persistence of pure T. oreades. Interestingly, 
when Crisp & Weston [8] excluded sclereids from their 
morphometric analysis, individuals sampled from Monga 
NP spanned the full range of morphological variation be- 
tween pure T. mongaensis and pure T. oreades. Their 
results were consistent with our genetic evidence but when 
viewed through the lens of a categorical distinction between 
plants possessing or lacking sclereids, morphological vari- 
ation at Monga NP was interpreted as circumscribing two 
sympatric, albeit sporadically hybridising species rather 
than a single, morphologically variable population [8]. 

Rather than a scenario of parapatric speciation, where 
hybrid zones continuously occur along adaptive clines at 
the distributional extremes of species [3], our findings are 
more supportive of a scenario similar to the one proposed 
for T. speciosissima and T. mongaensis. Temporal climatic 
fluctuations caused distributional changes that resulted in 
repeated cycles of secondary contact. As suitable habitat 
contracted, the T. oreades population at Monga NP 
remained isolated and was integrated by T. mongaensis. A 
history of allopatric speciation is unlikely to result in rapid 
and strong postzygotic reproductive isolation and conse- 
quently, complex hybridisation events can be expected in 
secondary contact zones [1]. However, further evidence for 
temporal changes in connectivity was needed in order to 
exclude a scenario of parapatric speciation. 

Bio-climatic evidence for temporal changes in 
connectivity 

A previous study on T. speciosissima suggested that the 
genetic differentiation measured among population groups 
representing coastal, upland and southern distributions 
was the consequence of cyclical distributional adjustments 
during the Quaternary [10]. An exploration of temporal 
changes in the availability of suitable habitat supports con- 
tractions into refugia during the LGM (Figure 2). In fact, 
Hesse et al. [64] suggested that local environmental condi- 
tions were so extreme during the LGM that upland areas 
were largely devoid of woody vegetation. 

Unexpectedly, whilst the temporal distributional variation 
that provides evidence for intra-specific allopatric differenti- 
ation in T. speciosissima is characterised by habitat contrac- 
tions during the LGM (also evident for its northern sister 
species, T. aspera), the circumstances are reversed for the 
two southern taxa (Figure 2). Current interglacial bio- 
climatic conditions confirm the existence of an inter-specific 
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distributional gap, while modelled LGM ENMs provide sup- 
port for expansions of suitable habitat. The LGM models 
support the scenario of a T. oreades x T. mongaensis hybrid 
population at Monga NP, currently disjunct from other 
T. oreades populations and likely to have originated from 
historical habitat expansions resulting in a now-vanished 
overlap zone. 

Overall, the bio-climatic models show that cyclical con- 
tractions/expansions of suitable habitat for all Telopea are 
likely to have generated temporally variable inter-specific 
distribution gaps consistent with allopatric differentiation. 
Contraction to small population size and the selection for 
even small adaptive differences would have enhanced the 
chance of faster stochastic allopatric differentiation [6]. 
Conversely, there is no evidence for the persistence of the 
continuous contact zones essential to a parapatric differen- 
tiation model [1]. 

Conclusions 

Although the traditional separation of speciation processes 
into allopatric and parapatric categories might not capture 
the full complexity of spatial relationships among taxa [65], 
our study shows that the evidence for temporal exclusion 
of gene flow can be found even outside obvious geograph- 
ical contexts (such as continental drift and island radiations 
for example). In Telopea, between-species differentiation 
has been regulated mostly by temporal changes in bio- 
climatic conditions that have repeatedly disrupted the con- 
tinuum in distribution of the various evolving lineages. 

We also show that it is possible to make significant pro- 
gress towards excluding parapatric speciation as a contrib- 
uting evolutionary process in a selected study system. 
Telopea includes allopatric taxa that are at different stages 
of speciation involving some ecological adaptation but no 
complete reproductive isolation. Niche differentiation does 
not imply the establishment of fixed reproductive barriers 
[3], and the combination of species-specific bio-climatic 
envelopes and hybrid zones suggests that selective filters 
shaped local adaptive distinctiveness as a secondary process 
rather than being the driver of speciation via the establish- 
ment of selection-based reproductive barriers. 

Current studies involving fine-scale coalescent analyses 
and transcriptome-based investigations of differential adap- 
tive potential among Telopea populations will bring further 
insights on speciation mechanisms operating in Telopea. 
Comparing genome-level patterns of diversity between 
populations can facilitate the identification of genomic 
regions that do not conform to the expectation of neutral 
demographic models. However, strong divergent selection 
at a small number of loci would not be expected in an allo- 
patric speciation scenario. Instead, analogous coalescent 
patterns across different loci could identify the signature 
of past periods of allopatry. For instance, targeting mul- 
tiple low recombination regions has been shown to be 



particularly useful in describing the contribution of allo- 
patric processes to sympatric divergences [66] . 
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